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^ ; ABSTRACT 

It is believed that turbulence may have a significant impact on star formation and the 
| dynamics and evolution of the molecular clouds in which this occurs. It is also known 

04 i that non-ideal magnetohydrodynamic effects influence the nature of this turbulence. 

We present the results of a numerical study of 4-fluid MHD turbulence in which the 
dynamics of electrons, ions, charged dust grains and neutrals and their interactions 
are followed. The parameters describing the fluid being simulated are based directly 
on observations of molecular clouds. We find that the velocity and magnetic field 
tIh ' power spectra are strongly influenced by multifluid effects on length-scales at least as 

large as 0.05 pc. The PDFs of the various species in the system are all found to be 
Q | close to log-normal, with charged species having a slightly less platykurtic (flattened) 

■ distribution than the neutrals. We find that the introduction of multifluid effects does 

not significantly alter the structure functions of the centroid velocity increment. 

c3 ; 
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\£) \ It is generally believed that molecular clouds are turbulent (see the reviews of Mac Low & Klessen 2004; Elmegreen & Scalo 

■ 2004). Observationally (e.g. Larson 1981; Brunt 2010) this turbulence appears to be highly supersonic with RMS Mach 
numbers of anything up to 20 or more. It is also thought that the Alfvenic Mach number is not significantly less than one, 
and may be much greater (e.g. Heyer & Brunt 2012). The amplitude of this turbulence makes it likely to be an important 
ingredient in both the dynamics of molecular clouds and the process of star formation (Elmegreen 1993; Klein et al. 2003) 

k> ■ and hence developing an understanding of this phenomenon is of considerable interest. 

Many authors have addressed the issue of MHD turbulence in the context of molecular clouds using both the ideal 
C$ • MHD approximation (Mac Low et al. 1998; Mac Low 1999; Ostriker et al. 2001; Vestuto et al. 2003; Gustaffson et al. 2006 
Glover & Mac Low 2007; Lemastcr & Stone 2008, 2009; Brunt, Federrath, & Price 2010a; Brunt, Federrath & Price 2010b 
Price, Federrath, & Brunt 2011) and, more recently, various flavours of non-ideal MHD (Oishi & Mac Low 2006; Li et al. 2008 
Kudoh & Basu 2008; Downes & O'Sullivan 2009, 2011). Three dimensional MHD turbulence involves the transfer of energy 
from an energy injection scale to ever smaller scales until the dissipation length scale of the system is reached. Given that in 
molecular clouds the lengthscales at which multifluid effects become important are much larger than the viscous lengthscale, 
it is inevitable that these effects will have an impact on the energy cascade. 

When ambipolar diffusion is included in MHD turbulence simulations it has been found that it causes greater temporal 
variability in the turbulence statistics (Oishi & Mac Low 2006; Li et al. 2008; Kudoh & Basu 2008). Although clearly of lesser 
significance (Wardle & Ng 1999; Downes & O'Sullivan 2009, 2011), the Hall effect is capable of inducing topological changes 
in the magnetic field which are quite distinct to any influence caused by ambipolar diffusion. In particular, it introduces a 
handedness into the flow which is of particular interest when considering the conversion of kinetic energy into magnetic energy. 
Researchers working on reconnection and the solar wind have studied the Hall effect in the context of turbulence and found 
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that, although the overall decay rate appears not to be affected, the usual coincidence of the magnetic and velocity fields seen 
in ideal MHD does not occur at small scales (Matthaeus et al. 2003; Mininni et al. 2006; Servidio et al. 2007). 

Downes & O'Sullivan (2009, hereafter Paper I) performed simulations of decaying, non-ideal MHD, molecular cloud 
turbulence incorporating parallel resistivity, the Hall effect and ambipolar diffusion. They found that the Hall effect has 
surprisingly little impact on the behaviour of the turbulence: it does not affect the energy decay at all and has very limited 
impact on the power spectra of any of the dynamical variables, with the exception of the magnetic field at very short length- 
scales. Downes & O'Sullivan (2011, hereafter Paper II) extended the results of Paper I to properly multifluid MHD turbulent 
decay. A system of three fluids, 1 neutral and two charged species, was simulated and, intriguingly, it appears that for the 
chosen system a full multifluid simulation is unnecessary and simple non-ideal MHD with temporally and spatially constant 
resistivities can give quite reliable results. 

In this paper we use the HYDRA code (O'Sullivan & Downes 2006, 2007) to investigate driven, isothermal, multifluid 
MHD turbulence in a system consisting of 4 fluids (1 neutral and 3 charged). The aim of this work is to determine the 
influence of multifluid effects on the behaviour of driven turbulence in molecular clouds. An added benefit of our properly 
multifluid approach is that we can self-consistently deduce the behaviour of the charged species and, in principle, make links 
with observations such as those of Li & Houde (2008). 

The structure of this paper is as follows: in section 2 we outline the equations and numerical method; in section 3 we 
specify our parameters, initial conditions and how we drive the turbulence in our simulations; section 5 contains the results 
of a resolution study demonstrating that the results we present are reasonably well converged; section 6 contains the results 
and analysis of our simulations and finally we draw conclusions from our work in section 7. We note here that we omit a 
discussion of the differences in line widths between the charged and the neutral species. The results associated with these 
differences require detailed discussion in their own right and are the subject of a forthcoming paper. 



2 THE MODEL 

Molecular clouds are, to a great extent, weakly ionised systems. This allows us to make several simplifying assumptions, as 
follows: 

• The bulk flow velocity is the neutral velocity 

• The majority of collisions experienced by each charged species occur with neutrals 

• The charged species' inertia is unimportant 

• The charged species' pressure gradient is unimportant 

Making these assumptions allows us to derive a relationship between the electric field and the current density. This generalised 
Ohm's law (e.g. Falle 2003; Ciolek & Roberge 2002) allows us to avoid having to calculate the electric field explicitly from 
the charge distribution. 

2.1 Equations 

Once we have our generalised Ohm's law, the system of equations for our weakly ionised flow is as follows (e.g. Falle 2003) 



31 Ul 



dt 



+ = 0, (1 < i < TV), (1) 

+ V • (puitti + a 2 p\) = J x B, (2) 



dB „ , „ „ , f (J ■ B)B J x B B x (J x B . , 

_+V.(«iB-Bui) = Vx r-o l B2 > + ri-g-+r 2 L_ L\ , (3) 



a t pi (e' + Ui x B) = -piPiKn{u! - Ui), 2 ^ i «£ N, (4) 
V-B = 0, (5) 
VxD = J, (6) 

AT 

^a lPl = 0. (7) 

i-2 

where i denotes the species (with i = 1 denoting the neutral species), pi, u%, a, B, I, J and E' are the mass density of species 
i, the velocity of species i, the isothermal sound speed, the magnetic field, the identity matrix, the current density and the 
electric field in the neutral rest frame, respectively. The parameters on and Kn are the charge-to-mass ratios and the collision 
coefficient with the neutrals for species i. 
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These equations describe conservation of mass for all species, the neutral momentum equation, the induction equation, 
force balance for the charged species, the absence of magnetic monopoles, Faraday's Law and charge neutrality. In the work 
presented here we set N = 4. 

The electric field in the frame of the fluid, E' , is calculated from the generalised Ohm's law for weakly ionised fluids (e.g. 

Falle 2003; O 'Sullivan & Downes 2006) and is given by 

E' = Eq + E\ + E2, (8) 
where 

E = (J -0,0)0,0, (9) 
Ei = J x 01, (10) 

E2 = —(J X 02) X 02, (11) 

using the definitions ao = foB, ai = fiB, 02 = /2-B, where fa = ^/ro/B, fi = ri/B, f2 = t/t^/B. ro, ri and T2 are the 
parallel, Hall and ambipolar resistivities respectively and are given by 

ro = —, (12) 

ri = whi> (13) 

- = (14) 

with the conductivities given by 
1 N 

Co = -gJ^^Pift, (15) 

i=2 
N 

^ = bI^ttwv (16) 

i=2 
AT 

1 \ - OLiPiPi , 

B^i+pr ( ' 

i=2 

where /3, is the Hall parameter for species i and is given by 

ft = (18) 



02 = 



2.2 Numerical method 

Equations (1) - (7) are a nonlinear system of PDEs and, in general, must be tackled using simulations. With modern techniques 
the system is relatively straightforward, with the exception of dealing with the right-hand side of equation (3). The method 
used by HYDRA to integrate these equations has been elucidated previously in the literature (e.g. O'Sullivan & Downes 2006, 
2007) so we only give a brief recapitulation of this here and refer the reader to these papers for a more in-depth discussion. 
The system is solved using an operator-split approach with three operators: 

• The equations governing the neutrals are solved using a standard second order, shock capturing Godunov scheme. The 
induction equation is also solved as part of this operator, but without incorporation of the non-ideal effects. The divergence 
of the magnetic field is controlled using the method of Dedner et al. (2002). 

• The magnetic diffusion terms are evaluated and applied using the super-time-stepping method (for the ambipolar and 
parallel resistivity terms) and the Hall Diffusion Scheme (for the Hall term). 

• The velocities of the charged species are determined from equation (4) and the distribution of their densities are updated 
using a second order upwind scheme. 

The whole system of operators is stepped forward in time using Strang operator splitting in order to ensure second order 
accuracy in time. 

Since molecular cloud material has very low kinematic viscosity (e.g. Jones & Downes 2012), flows, these systems have a 
high Prandtl number and thus we can consider them to be effectively inviscid. As with all higher-order schemes, our simulations 
incorporate artificial viscosity which smears discontinuities over around 3-5 grid zones. This artificial viscosity length scale is 
much less than that over which the multifluid effects act (see, for example, figure 2) and thus our multifluid effects dominate 
at all scales above around 3-5 grid zones. Formally, the Prandtl number of the highest resolution simulations in this work 
is around 15. 
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Name 


Resolution (Ng) 


Physics 


mf-64 


64 


Multifluid MHD 


mf-128 


128 


Multifluid MHD 


mf-256 


256 


Multifluid MHD 


mhd-256 


256 


Ideal MHD 



Table 1. Definition of the nomenclature used to refer to the various simulations presented in this work. 

The method is explicit and scales extremely well on massively parallel systems, displaying strong scaling from 8192 up 
to 294,912 cores with 73% efficiency on the JUGENE BlueGene/P system in Jiilich, Germany. 

3 NUMERICAL SET-UP 

We run our turbulence simulation in a cube of side 0.2 pc with periodic boundary conditions on all faces. The density is initially 
uniform and has the value of 10 4 cm~ 3 , assuming that the average neutral particle has a mass of 2.33 m p (where m p is the 
mass of the proton). The magnetic field is also initially uniform of strength 20 /iG respectively. The (isothermal) sound-speed 
is chosen to be 1.88 x 10 4 cms -1 or, equivalently, the temperature is 10 K. The four fluids in the system are initially at rest. 
These parameters were chosen to be representative of molecular clouds on scales of 0.2 pc and are based on the data in Table 
1 of Crutcher (1999). The grid used is uniform and made up of N x X N y x N z grid points with N x = N y = N z = iV g . We 
emphasise here that this work is directed at understanding the nature of multifluid MHD turbulence. A physical system with 
the properties just described will have a mass just below the Jeans mass and so may well be subject to dynamically important 
self-gravity. However, it seems reasonable that the optimal way of approaching understanding the behaviour of molecular 
clouds is to first attempt to understand the turbulence, and then to introduce self-gravity. Introducing all the physical effects 
at once runs the risk of producing results which, while possibly valid, are very difficult to interpret. 

The properties of the charged species are based on the data in Wardle & Ng (1999). We make the approximation that 
the dust density is 1% that of the neutrals. The metal ions are assumed to be singly ionised and to have an average mass of 
24 m p and a mass density of 1.27 x 10~ 25 gcm~ 3 . The electron density is then chosen to ensure local charge neutrality. This 
set-up results in an ionisation fraction of approximately 3 x 10~ 7 and the system is therefore clearly weakly ionised. 

In order to drive turbulence we add velocity increments to the neutral velocity at each time-step. The incremental velocity 
field, Su, is defined in a manner similar to Lemaster & Stone (2009). Each component of Su is generated from a set of waves 
with wave numbers, k = |fe|, 3 ^ k ^ 4. The amplitudes of these waves are drawn from a Gaussian random distribution with 
mean 1.0 and deviation 0.33 while the phases of the waves are drawn from a uniform distribution between and 2-7T. The waves 
are chosen such that the forcing is solenoidal. It is known that the precise nature of the forcing will influence the turbulence 
statistics (e.g. Fcderrath et al. 2011). In this paper we focus on gaining some insight into the influences of multifluid effects 
on turbulence and leave a study of multifluid MHD turbulence under differing driving to a later date. 

Finally, to ensure a constant energy injection rate and zero net momentum injection, a quadratic equation must be solved 
to normalise Su. The resulting normalised incremental velocity field is then added on to the neutral velocity field. The rate of 
injection of energy is the same for all simulations in this work and is given by po ^i a ii = 200 where po is the initial (uniform) 
density, L is the length of side of the computational domain and a is the (isothermal) sound speed. In our simulations this 
energy injection rate yields steady state turbulence with an RMS Mach number of around 4.5 and an Alfvenic Mach number 
of between 3 and 3.5. 

Table 1 contains the definition of the nomenclature of the various simulations used in this work. 



4 ANALYSIS 

There are two different, but related, aspects to analysing the multifluid simulations presented here. The first is associated 
with understanding the dynamics of the bulk flow, while the second is focused on the interplay, and differences, between the 
behaviour of the various species in the simulation. In both cases, however, we can use analysis techniques which are broadly 
similar. 

4.1 Power spectra 

One of the most popular things to do when faced with a set of results from turbulence simulations is to take power spectra. 
These spectra give information on the relative size of structures formed by turbulence in the system. In incompressible 
hydrodynamic turbulence this is easily visualised as, for example, the energy residing in vortices of differing sizes. Standard 
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Simulation 


M B 


M A 


Velocity power spectrum 


mf-64 


4.38 


2.77 


-2.54 ±0.26 


mf-128 


4.47 


3.01 


-2.40 ±0.23 


mf-256 


4.48 


3.07 


-2.40 ±0.05 


mhd-256 


4.41 


3.17 


-1.81 ±0.05 



Table 2. The time-averaged RMS sonic Mach and Alfven numbers and the slope of the velocity power spectra for the simulations 
presented in this work. The velocity power spectrum is fitted over the range [4 : N g /16] in each case with the exception of mf-64 where 
the range is [4:8]. N g is the number of grid zones along a side of our cube. The error quoted is the formal error from the least squares 
fitting. 



Kolmogorov theory predicts that, in the absence of driving or dissipation, the power spectrum will take on a power-law 
form, and the exponent of this power law can be predicted in a relatively straightforward manner. For incompressible MHD 
turbulence the situation is a little more complicated, while for compressible MHD turbulence it becomes more complicated 
again. However, the power spectra of the fluid variables still provide useful physical insight into the type of turbulence 
occurring. 

The power spectra presented here are calculated by taking a 3D Fourier transform of the relevant fluid variable. The total 
power residing in a spherical shell of inner radius k and outer radius k + dk is then calculated and it is this quantity which we 
present in our power spectra. When discussing wave numbers we use the notation k = ^k where k is the usual wave number. 
Then k is the number of waves which fit into the computational domain for that value of k. In calculating the power spectra 
we use dk = 1. The power spectra are taken for different times and the results presented here are the time-averaged spectra. 
For all simulations the spectra are averaged over t — [0.2, 0.52] sound crossing times, or [0.9, 2.34] flow crossing times. 



4.2 The centroid velocity 

The statistics of the centroid velocity are, in principle at least, an observable quantity and are therefore of some physical 
interest (e.g. Hily-Blant, Falgarone, & Pety 2008; Federrath et al. 2010). Suppose we orient our simulation cube such that 
the z direction is along the line of sight, leaving the x and y directions in the plane of the sky. The centroid velocity is then 
defined as 

c(x,v) = - — n 73 ( 19 ) 

The density weighting used here to calculate the average line-of-sight velocity can be thought of as some approximation to 
what happens when observations are made and the average line-of-sight velocity is actually the emission-weighted average. 
Having defined our centroid velocity we can then consider the centroid velocity increment, defined by 

5G(t) = {C(t)-C(t + 1)) (20) 

where r = (x,y) and 1 is a vector of length I and the average is calculated over all possible directions of 1 in the (x,y) plane. 
Finally we introduce the associated centroid velocity increment structure functions, defined by 

CVISF p (0 =< |5C ; (r)| p > r (21) 

where p is the order of the structure function. These quantities are discussed in Sect. 4.2. 



5 RESOLUTION STUDY 

We performed identical simulations with differing resolutions in order to determine what aspects of the results from our highest 
resolution simulation are well converged. The resolutions used were N s = 64, 128 and 256 (see Sect. 3). Table 2 contains a 
summary of our results from the resolution study. 

Figure 1 contains plots of the RMS Mach number as a function of time for each of the simulations in our resolution 
study. It is clear that after an initial period of time during which the turbulence is building up, each simulation reaches a 
quasi-steady state with an RMS turbulent Mach number of between 4.3 and 4.5. 

The highest Mach number reached occurs, perhaps unsurprisingly, in the 256 simulation. However, after this initial 
maximum, the average Mach number in each case is rather similar, bearing in mind the stochastic nature of the simulations, 
with the 128 3 and the 256 3 simulations giving results within about 0.2% of each other. Simulation mf-64 does have a somewhat 
lower RMS Mach number (see Table 2), indicating that at this low resolution the global energy dissipation rate is significantly 
influenced by numerical viscosity. However, the similarity of the Mach numbers for mf-128 and mf-256 indicate that the global 
energy dissipation rate is well captured with a resolution of 256 3 . 
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Figure 1. Plot of the RMS Mach number as a function of time for each of the simulations in the resolution study. 



U 
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k 

Figure 2. Plot of the time-averaged, compensated neutral velocity power spectrum for each of the simulations in the resolution study. 
See also Table 2. The spectra are compensated by fc 2,2 . 

Another measure of whether or not turbulence simulations have converged is the slope of the power spectrum of the 
velocity. Figure 2 contains plots of the compensated power spectrum of the neutral velocity for each of the three simulations 
in the resolution study. While we generally expect this power spectrum to behave as a power law, there will be some value 
of k at which the spectrum will steepen due to numerical dissipation effects. In order to get a picture of the range of k over 
which we can trust our results, we have fitted the power spectrum with a power-law over various different ranges. The driving 
occurs at 3 ^ k ^ 4 and we find that between k = 4 and k « (N g being the number of grid zones along a side of our cube) 
our slopes are very similar, keeping in mind the formal errors, across all resolutions (see Table 2). These results give support 
to the idea that these simulations have resolved the turbulent cascade, at least up to a value of k of around -j^ zones. Even 
so, when interpreting the results in a quantitative fashion, the size of the formal errors must be kept in mind and the best 
way to reduce these is through using higher resolution simulations. The results of Federrath et al. (2010) suggest that caution 
must be used when interpreting power spectra at high k. In their particular work where FLASH3 was used they found the 
power spectra to be susceptible to numerical dissipation below length scales of about 30 grid zones. Our resolution study here 
demonstrates that the equivalent length scale, for velocity power spectra, for HYDRA is about 10 - 15 grid zones. 



6 RESULTS 

Following the results of our resolution study we are reasonably confident that our N g = 256 simulation is well enough resolved 
to draw some useful physical conclusions about the behaviour of driven, multifluid MHD turbulence. 

Figure 3 contains plots of the (neutral) column density for the multifluid and ideal MHD simulations. It is apparent that 
there is somewhat more small scale structure in mhd-256 than in mf-256. The maxima and minima of the column density are 
broadly similar in each case. 
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5 - 
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Figure 3. Distribution of the column density for mf-256 (upper panel) and mhd-256 (lower panel) at time t = 0.5 t c , assuming the line 
of sight is along the z axis. The units on the x and y axes are grid zones and the density is in units of 10 4 cm -3 . 



6.1 Energy dissipation 

Figure 4 contains plots of the behaviour of the mean, mass-weighted RMS Mach and Alfven numbers (see Lemaster & Stone 
2009; Downes & O'Sullivan 2011). It is clear that the Mach number reached in each case is very similar. Table 2 demonstrates 
this in quantitative terms with the difference in time-averaged RMS Mach number after t = 0.2 i c being around than 1.6%. 
Therefore, at these length scales and with this numerical resolution, we do not detect any difference in the energy dissipation 
rate due to multifluid effects. Note that the RMS Mach number is actually a simple measure of the RMS velocity on the 
computational grid as the simulations are isothermal. On the other hand, the RMS Alfven number is a more complicated 
quantity, given that the Alfven speed is a spatially and temporally varying quantity. It is clear that the difference in the 
Alfven numbers between mhd-256 and mf-256 are rather small, with a difference of only around 3%. This is quite interesting, 
given that ambipolar diffusion does diffuse magnetic energy and hence would be expected to lower the Alfven speed generally. 
Hence, taken together, the similarities between mf-256 and mhd-256 with respect to their Mach and Alfven numbers indicates 
that the Alfven speed is not, in fact, lowered significantly on a global scale by multifluid effects. 

The results presented in Papers I and II imply that the energy dissipation rate is higher in multifluid MHD systems. We 
would expect, then, that the mean Mach number would be lower in mf-256 than in mhd-256. However, it is clear from Table 
2 that this is not the case. If, instead of the time-averaged mass-weighted RMS Mach number, we calculate the time-averaged 
RMS Mach number without mass- weighting, we find that the RMS Mach numbers are 4.18 and 4.10 for mhd-256 and mf-256 
respectively. This is a little more in line with what we would naively expect although, on the face of it, it still does not indicate 
significantly increased energy dissipation due to multifluid effects. When considering this issue, it is significant that the main 
contribution to the RMS Mach numbers in each of mhd-256 and mf-256 comes from large scale motions - i.e. motions close to 
the driving scale. Since the multifluid effects operate at relatively short length scales it, perhaps, is not too surprising that the 
RMS Mach number is not effected much by the inclusion of these effects. When reconciling this result with those of Papers 
I and II it is worth recalling that these latter results were calculated for decaying turbulence. As turbulence decays energy 
initially injected at large scales cascades to ever smaller scales. Thus decaying multifluid turbulence will exhibit a faster Mach 
number decay as, in contrast to driven turbulence, the main contribution to the measured Mach number will come from ever 
smaller scales, and thus be more influenced by multifluid effects. 

There is a further effect to take into account. Closer examination of the data from mf-256 and mhd-256 reveals that, while 
the kinetic energies in each simulation are very similar, the magnetic energy in the multifluid simulation is much decreased (see 
Figure 5). This reduction in the magnetic energy is due to the presence of ambipolar diffusion and shows that the multifluid 
effects do have a significant impact on the energetics of the system. The fact that the RMS Mach number is not strongly 
influenced by multifluid effects is partly a reflection of the fact that the energy in the system is dominated by kinetic energy. 
Some of this kinetic energy is converted into magnetic energy but, even in mhd-256, this is a relatively small proportion: on 
average the magnetic energy is around 40% of the kinetic energy. The fact that this converted energy is then lost faster in 
mf-256 has relatively little impact on the large reservoir of kinetic energy, and therefore has a minor influence on the RMS 
Mach number. On a side note, it is interesting to note the results of Fcdcrrath et al. (2011) in which studies of dynamo action 
in MHD simulations of turbulence. These authors found that, as might be expected, solenoidally forced turbulence induces 
stronger dynamo action. The relevant saturation levels for the ratio of magnetic to kinetic energy from that study are around 
0.03, considerably smaller than the results presented here. However, the significant differences in terms of initial conditions 
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mf-256: Mach number 
mhd-256: Mach number 
mf-256: Alfven Number 
mhd-256: Alfven Number 




Time (crossing times) 

Figure 4. Plot of the RMS Mach and Alfven numbers as functions of time for mf-256 and mhd-256. 




Time (sound crossing times) 

Figure 5. Plot of the total perturbed magnetic energy in the computational domain as a function of time for mf-256 and mhd-256. 

between this latter study and the one presented here make it difficult to draw any firm conclusions about the origins of the 
differences between the results. 

In the case of decaying turbulence, the magnetic energy becomes approximately equal to the kinetic energy as the decay 
process continues (see Figure 6 of Paper I) and hence the multifluid effects can have a more significant impact on the overall 
energetics of the system, resulting in a faster decay rate. 

6.2 Power spectra 

We now turn to the power spectra of the various quantities of interest in these simulations. The nature of these spectra will 
give us some idea of the influence of the multifluid effects on the turbulence. As noted previously, we time-average all power 
spectra discussed here over the time-interval [0.2 t c : 0.52 t c ]. Table 3 contains the exponents of the power spectra, fitted 
assuming the spectrum goes as k a . The fitting is carried out over the same range as that in Table 2 and so we are reasonably 
confident that the exponents we find are not significantly influenced by numerical resolution effects. 



Table 3. Exponents of the power spectra of the various variables in the mf-256 and mhd-256 simulations. The spectra were fitted as 
power laws in the range 4 ^ k < 





mhd-256 


Neutrals 


mf-256 
Electrons 


Dust 


Ions 


Density 


-0.60 


-0.83 


-0.99 


-1.11 


-0.99 


Velocity 


-1.81 


-2.40 


-2.17 


-2.19 


-2.17 


B-field 


-1.73 




-2.26 
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mhd-256 
mf-256 




Figure 6. Time-averaged compensated, spherically integrated velocity power spectra for the neutrals in mf-256 and the bulk flow in 
mhd-256. The spectra are compensated with fc 18 . 




Neutrals 
Electrons 
Dust 



Figure 7. Time-averaged, compensated, spherically integrated velocity power spectra for the neutrals, electrons and dust in mf-256. 
The spectrum for ions is virtually identical to that for electrons and is omitted. 



6.2.1 Velocity power spectra 

The compensated velocity power spectra for the neutrals in mf-256 and the bulk fluid in mhd-256 are shown in figure 6. The 
power is compensated with fc 1,8 which is the measured exponent of the velocity spectrum of mhd-256 (see Table 3). It is clear, 
as indicated also in Table 3, that the multifluid effects modify the velocity power spectra of the neutrals significantly, even on 
rather large length scales. 

Figure 7 contains plots of the spherically integrated, compensated velocity power spectra for the neutrals, electrons and 
dust in mf-256. We have omitted a plot of the spectrum for the ions as it is virtually identical to that for the electrons. It is 
clear that the velocity power spectra for the neutrals and charged species do differ to some extent on these length scales (i.e. 
all scales up to 0.05 pc). The spectrum for the charged species is somewhat harder than the neutrals. This is in line with what 
would be expected in a multifluid system affected by ambipolar diffusion: as ambipolar diffusion becomes more important the 
neutrals behave more like a hydrodynamic system which characteristically has a softer velocity power spectrum than an MHD 
system. The charged species remain more affected by the magnetic field and therefore retain more strongly the characteristics 
of MHD turbulence. Another way of looking at this is to observe that, from the data in Table 3, the power spectra of the 
velocity for the charged species have exponents which lie between that for the neutrals and that for ideal MHD. Again, this 
is what would be expected given the more direct link between the behaviour of the charged species and that of the magnetic 
field. At very high k there is a slight difference between the dust and electron power spectrum, almost certainly due to the 
Hall term in the induction equation. However, we should refrain from attaching quantitative physical significance to this as 
the length scales involved at these values of k are too short to be reliably resolved in the simulation. 
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Figure 8. Plots of the time-averaged power spectra for the neutral density in mf-256 and the density in mhd-256. The spectra are 
compensated using k 06 (see Table 3). 
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Figure 9. Plots of the the time-averaged, normalised power spectra of the density for all species in mf-256. The spectra are normalised 
and compensated by k 06 for ease of comparison (see Table 3). The power spectrum for the ions is not plotted as it is identical to that 
for electrons. 



6.2.2 Density power spectra 

Density power spectra, while not as dynamically relevant in supersonic turbulence as velocity power spectra, nonetheless 
provide an insight into the size of structures formed by the turbulence (e.g. Ossenkopf & Mac Low 2002). This can also be 
investigated in a different way using the probability density function (PDF) of the density which we discuss later (Sect. 6.3). 
Figure 8 contains plots of the time-averaged (neutral) density power spectra for mf-256 and mhd-256. 

There is less power in the bulk density of mf-256 than mhd-256 for all k ^ 10 as one might expect, given the similarity 
between the multifluid effects and viscosity and also given the behaviour of the velocity power spectra. Interestingly, there is 
marginally less power in mhd-256 than mf-256 at k < 10. This difference is rather small and is likely to be due to diffusive 
effects causing an inverse cascade, as would be expected since they smooth out small scale structures into larger scale ones. 

While there is somewhat less power in the density variations in mf-256 than mhd-256, there are striking differences 
between the power in the density variations of the neutral and charged species in mf-256 (see Figure 9). Hence we have a 
situation where the neutral power spectrum is more qualitatively similar to that obtained from ideal MHD, while the charged 
species (which are more directly influenced by the magnetic field) have significantly less power at virtually all length scales. A 
careful comparison of figure 10 and the upper panel of figure 3 reveals that, even by eye, there is less structure in, for example, 
the dust distribution than the neutral distribution. 

The inclusion of multifluid effects leads to the introduction of terms in the induction equation which are diffusive so it 
is not a surprise that there is less structure in a multifluid MHD system than an ideal MHD one. It is also, perhaps, not 
surprising that this lack of structure is particularly marked in the charged species which are more directly affected by the 
dynamics of the magnetic field. 
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Figure 10. Distribution of the column density for the dust in our mf-256 simulation at time t = 0.5 t c , assuming the line of sight is along 
the z axis. The units on the x and y axes are grid zones and the density is in units normalised by the initial neutral density. Comparison 
with the upper panel in figure 3 indicates there is less small-scale structure in the dust distribution. 
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Figure 11. Time-averaged, compensated, spherically integrated magnetic field power spectra for mf-256 and mhd-256. Each power 
spectrum is compensated with fc 7 / 4 . 



6.2.3 Magnetic field power spectra 

We now turn our attention to the power spectra of the magnetic fields in each simulation (see Fig. 11). In keeping with 
the velocity and density power spectra (Figs. 8 and 6), there is less structure on all scales, even at the driving scale, in the 
multifluid system than in the ideal MHD system. This is due to the presence of ambipolar diffusion which, being a standard 
diffusion term in the induction equation, will smear out any small scale structure. This result implies that multifluid effects 
can have an impact on the nature of molecular cloud turbulence even at fairly large length scales - much larger than those 
inferred from, for example, the observations of Li et al. (2008, 2010). We address this apparent contradiction in a forthcoming 
paper. 

It is worth noting again that the Hall effect is present in this system, albeit at a much lower level than ambipolar diffusion. 
As in Paper I and Paper II we find little evidence of the Hall effect at length scales which are numerically resolved in our power 
spectra. At very short length-scales (k ^ 100) we can see an upturn in the power spectrum for the magnetic field and the 
charged species (Figure 7) and this is almost certainly the influence of the Hall effect where, as expected, we see a decoupling 
at short length scales between the neutrals and the charged species, as well as the magnetic field. The phenomenon of the 
Hall effect having somewhat less impact than would be expected was also discussed in Jones & Downes (2012) who deduced 
that the presence of ambipolar diffusion can effectively short circuit the Hall current, thereby inhibiting the Hall effect even 
in systems with significant Hall resistivity. 
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Figure 12. Time-averaged PDF of the normalised bulk density for mf-256 and mhd-256. 

Table 4. Statistics of the probability distribution functions of the logarithm (to the base 10) of the various densities in the simulations 
presented here. 





mhd-256 




mf-256 










Neutrals 


Electrons 


Dust 


Ions 


Standard deviation 


0.46 


0.44 


0.36 


0.34 


0.36 


Skewness 


-0.050 


-0.12 


0.15 


0.082 


0.15 


Kurtosis 


-0.11 


-0.071 


-0.005 


0.057 


-0.005 



6.3 Probability Density Functions 

Probability density functions (PDFs), in particular of the mass density, are particularly interesting when considering the 
influence molecular cloud turbulence may have on star formation. It has long been known that isothermal, supersonic turbu- 
lence results in a log-normal PDF of the density (e.g. Vazquez-Semadeni 1994; Passot & Vazquez-Semadeni 1998). This is due 
simply to the fact that the density of a particular fluid element is a product of its initial value and all of the enhancements 
at each of the shocks it has encountered. Since the shock strengths are randomly distributed, the density is thus the result 
of a product of random variables and hence its distribution is log-normal. Several researchers have shown that the addition 
of gravity to the physics of the system adds a long, high mass tail to the PDF (e.g. Klessen 2000; Dib & Burkert 2005), 
presumably due to the influence of gravity which will produce accumulations of mass in small regions which would not be 
expected to happen to the same extent with isothermal turbulence. This would lead to a positive skewness of the density 
distribution or, to say it another way, the distribution being skewed right. 

Figure 12 contains plots of the PDF of the logarithm of the density of the bulk flow of mhd-256 and mf-256. Table 4 
details the standard deviation, skewness and excess kurtosis of the distributions of all of the densities in our simulations. 
Clearly all of the densities in each of mf-256 and mhd-256 have very close to a log-normal distribution. The logarithm of the 
bulk mass density is very slightly platykurtic in both mf-256 and mhd-256, while the distributions of the charged densities 
are a mix of leptokurtic and platykurtic. However, these differences in the kurtosis are rather minor. The skewness of these 
distributions is also small, indicating again how close to log-normal the PDF of the density of the bulk flow is. 

Examining Figure 13 and Table 4 we see that the charged species are slightly less platykurtic (flattened) than the neutrals 
in mf-256. We can understand this as follows: the magnetic field has, of course, a pressure associated with it which will resist 
both compression and rarefaction. Since the charged species are, by definition, more closely tied to the magnetic field than the 
neutrals we expect that the neutrals will undergo more severe compressions and rarefactions than the charged species. Thus 
the distribution of the charged species densities' should be more peaked around their mean than the neutrals. The standard 
deviation of our neutral species implies a value for the constant of proportionality between the Mach number and standard 
deviation (the "b" parameter in, for example, Fcderrath, Klessen, & Schmidt 2008) of 0.31, in keeping with the solenoidal 
forcing used in the work presented here. That the charged species have a higher standard deviation is unlikely to be of direct 
importance to models of the Initial Mass Function as the charged species are insignificant in terms of their mass, and hence 
their direct contribution to gravitational collapse. That the charged species are less abundant in high and low density regions 
than in intermediate density regions does have implications for the behaviour of the resistivities in these regions and hence for 
the influence of ambipolar diffusion and fragmentation. A detailed study of this is, however, beyond the scope of this paper. 
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Figure 13. Time-averaged PDF of the normalised densities of the neutrals, electrons and dust in mf-256. The ion density PDF is not 
shown as it is identical to that of the electrons. 





Figure 14. Plots of the first order (top panel) and second order (bottom panel) structure functions of the centroid velocity increment 
for the neutral species in mf-256 and the bulk flow in mhd-256. Both these clearly follow power laws in I up to lengths of around one 
quarter of the computational domain. 



6.4 Centroid velocity statistics 

We now turn to the statistics associated with the centroid velocity (see Sect. 4.2). Figure 14 contains plots of the first 2 
structure functions of the centroid velocity (equation 21). Table 5 contains the absolute and relative exponents of these 
structure functions found by the fit CVISF p (i) oc fa. From both Figure 14 and Table 5 it is clear that the introduction of 
multifluid effects appears not to influence the centroid velocity structure functions of the neutral flow significantly. 

It is worth noting that the absolute exponents of the structure functions differ somewhat from those given in Fedcrrath et al. 
(2010) with a difference of around 15% for the p — 1 structure function. We performed a further ideal MHD simulation with 
a resolution of 512 3 to assess whether resolution issues may be affecting our results in this regard and we were able to confirm 
our results for the exponents. Since the simulations of Fedcrrath et al. (2010) were hydrodynamic we conclude that the reason 
for the differences in our results is due to the presence of dynamically significant magnetic fields. 
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Simulation 


mf-256 


mhd-256 


Ci 


0.32 (0.71) 


0.32 (0.69) 


C2 


0.65 (1.43) 


0.65 (1.39) 


C3 


1.00 (2.19) 


1.00 (2.14) 



Table 5. The exponents of the CVISF P (2), normalised to p = 3 for p = 1, 2, 3. The numbers in brackets are the absolute exponents. 
7 CONCLUSIONS 

We have presented the results of 4-fiuid, MHD turbulence in molecular clouds. The parameters chosen were matched to the 
observations of molecular clouds reported in Crutcher (1999). We performed a resolution study to ensure that the results 
presented, in terms of both the power spectra and the RMS Mach number, were not strongly influenced by resolution effects. 

We found that the RMS mass-weighted Mach number was not significantly decreased with the inclusion of multifluid 
effects, in spite of multifluid effects providing an effective pathway for removal of magnetic energy. The overall energy dissi- 
pation, then, appears not to be significantly affected by multifluid effects. This phenomenon is a result of the fact that the 
majority of the energy in a driven, turbulent system with the properties of a molecular cloud resides in kinetic energy and so 
an enhanced removal of magnetic energy does not significantly effect the global energetics of the system. 

We found that both the velocity power spectra of the neutrals and the magnetic field power spectra were strongly 
influenced by multifluid effects at all length-scales up to the driving scale (0.05 pc). In multifluid MHD, the velocity power 
spectrum of the neutrals is found to have considerably less structure at all length-scales than those of the charged species. 
The electron and ion fluids behave almost identically to each other, while the dust displays some slight differences with the 
other charged species. This latter point arises largely from the requirement of local charge neutrality and the fact that the 
electrons and ions are the dominant charge carriers in the system. We also find, at very short length-scales, some signs of the 
Hall effect in the power spectra of the charged species. 

Thus we find, in common with the results presented in Papers I and II that multifluid effects are important up to quite 
large length-scales. This result appears to be in conflict with the observations of Li & Houde (2008) and we deal with this 
apparent contradiction in a forthcoming paper. 

The differences in the power spectra of the bulk densities in the multifluid and ideal MHD systems are not as dramatic as 
those in the other variables, although there is somewhat less power in the distribution of the neutral density in the multifluid 
simulation than the bulk density in the ideal MHD one. Interestingly, there is less power in the spectrum of the charged 
species density than the neutral species in the multifluid system. This is the opposite of what is found for the velocity power 
spectra. 

For all species, and for all cases, the density PDFs are found to be very close to log-normal. The neutral species in the 
multifluid simulation has a slightly less platykurtic (flattened) distribution than that in the ideal MHD case. The charged 
densities in the multifluid simulation have, in turn, less platykurtic distributions than the neutrals in that simulation with 
the dust actually having a leptokurtic distribution. Each of these results is explained in terms of the ambipolar diffusion 
experienced by the multifluid system, and the fact that the charged species are more directly influenced by the magnetic field 
than the neutrals. 

Finally, we find that multifluid effects do not significantly influence the centroid velocity increment structure functions. 
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